begin
  
  fi = addfile("../uv_regional4_mcv.nc", "r")
  power_div  = fi->power_div
  power_vor  = fi->power_vor
  wavelength = fi->wavelength
  wavenumber = fi->wavenumber 

  data = new((/3, dimsizes(power_div)/), double)
  data(0,:) = power_div + power_vor
  data(1,:) = power_vor
  data(2,:) = power_div
  
  
 ; Lindborg(1999) formula 
  d1 = 9.1e-04
  d2 = 3.0e-10
  E_lindborg = new((/8*dimsizes(wavenumber)/), double)
  wavenumber2 = fspan(1.0e-06, 1.0e-02, 8*dimsizes(wavenumber))
  pi = get_pi("double")
  wavelength2 = 2 * pi / wavenumber2 * 1.e-03 ; km
  E_lindborg = d1 * wavenumber2^(-5.0/3.0) + d2 * wavenumber2^(-3.0)
 ; grid spacing
  radius = 6.37122e03 ; km
  resolution = 0.25
  dx = 2 * pi * radius / 360 * resolution

  wks = gsn_open_wks("ps", "../power_uv_mcv4")
  res1                    = True
  res1@gsnFrame           = False
  res1@trXLog             = True
  res1@trYLog             = True
  res1@trXMinF            = 1.0e-06 ;min(wavenumber)
  res1@trXMaxF            = 1.0e-02 ;max(wavenumber)
  res1@trYMaxF            = 1.e09
  res1@trYMinF            = 1.e00
  res1@xyLineThicknesses  = (/4.0, 4.0, 4.0/)
  res1@tiXAxisString      = "Wavenumber (radians m~S~-1~N~)"
  res1@tiYAxisString      = "E(k) (m~S~3~N~/s~S~2~N~)"
  res1@xyDashPatterns     = (/0, 0, 0/)
 ; res1@xyLineColors       = (/"black", "blue", "red"/)
  res1@xyLineColor        = "black"

  res2                = True
  res2@gsnFrame       = False
  res2@trXLog         = True
  res2@trYLog         = True
  res2@trXMinF        = min(wavelength2)
  res2@trXMaxF        = max(wavelength2)
;  res2@trYMinF        = res1@trYMinF
;  res2@trYMaxF        = res1@trYMaxF
  res2@trXReverse     = True
  res2@tmXTMode       = "Explicit"
  res2@tiXAxisString  = "Wavelength (km)"
  res2@tmXUseBottom   = False
  res2@tmXTLabelsOn   = True
  res2@xyDashPatterns = res1@xyDashPatterns
;  res2@xyLineColors   = res1@xyLineColors 
;  res2@gsnXRefLine    = (/2 * dx, 4 * dx/)
;  res2@gsnXRefLineDashPatterns = (/0, 1/)
;  res2@gsnXRefLineThicknesses = (/4, 4/)
;  res2@gsnXRefLineColors = (/"cyan", "cyan"/)
;  gsn_text_ndc(wks, "2~F33~D~F~", 0.5, 0.5, res2)

  plot = gsn_csm_x2y(wks, wavenumber2, wavelength2, E_lindborg, res1, res2)
  res1@xyLineColor    = "red"
  res1@tmXTOn         = False
  plot = gsn_csm_xy(wks, wavenumber, power_div+power_vor, res1)
  
  ;add legend
  lgres                     = True
  lgres@lgLabelFontHeightF  = 0.1
  lgres@vpWidthF            = 0.25
  lgres@vpHeightF           = 0.1
  lgres@lgPerimOn           = True
  lgres@lgItemOrder         = (/0,1/)
  lgres@lgDashIndexes       = (/0,0/)
  lgres@lgLineColors        = (/"black", "red"/)
  lgres@lgLineThicknessF    = 4.0
  labels                    = (/"Lindborg (1999), eqn 71", "GRAPES"/)
  legend                    = gsn_create_legend(wks, 2, labels, lgres)
  amres                     = True
  amres@amZone              = 0
  amres@amParallelPosF      = -0.25 ; left-right
  amres@amOrthogonalPosF    = 0.4 ; up-down
  annota = gsn_add_annotation(plot, legend, amres)
  
  draw(plot)
  frame(wks)
end

